from util import *
from glob import glob
/usr/local/lib/python3.10/dist-packages/geopandas/_compat.py:123: UserWarning: The Shapely GEOS version (3.11.1-CAPI-1.17.1) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow. warnings.warn( /home/nyou045/git/CliffErosion_Projections/util.py:3: UserWarning: Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas will still use PyGEOS by default for now. To force to use and test Shapely 2.0, you have to set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas: import os os.environ['USE_PYGEOS'] = '0' import geopandas In a future release, GeoPandas will switch to using Shapely by default. If you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html). import geopandas as gpd # Geospatial data
df = load_AOIs()
df
| Taranaki | AOI | SSP 4.5 (p50) | SSP 4.5 (p83) | SSP 8.5 (p50) | SSP 8.5 (p83) | Rate SSP 4.5 (p50) | Rate SSP 4.5 (p83) | Rate SSP 8.5 (p50) | Rate SSP 8.5 (p83) | match | match_score | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 7 | NORTH | TongaporutuRiver | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | TongapurutuRiverCliffs | 93.750000 |
| 11 | SOUTH | HangatahuaRiver_South | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | HangatahuRiver_South | 97.435897 |
| 21 | SOUTH | Rahotu | 0.58 | 0.78 | 0.84 | 1.10 | 0.0058 | 0.0078 | 0.0084 | 0.0110 | Rahotu | 100.000000 |
| 20 | SOUTH | Pihama | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Pihama | 100.000000 |
| 19 | SOUTH | OpunakeBeach | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | OpunakeBeachCliffs | 100.000000 |
| 18 | SOUTH | OhaweBeach | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | OhaweBeach | 100.000000 |
| 17 | SOUTH | Oeo | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Oeo | 100.000000 |
| 16 | SOUTH | Manutahi | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | Manutahi | 100.000000 |
| 15 | SOUTH | ManaBay | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | ManaBayCliffs | 100.000000 |
| 14 | SOUTH | KaupokonuiBeach | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | KaupokonuiBeach | 100.000000 |
| 13 | SOUTH | Kakaramea | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | Kakaramea | 100.000000 |
| 12 | SOUTH | Hawera_WaihiBeach | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | Hawera_WaihiBeach | 100.000000 |
| 0 | NORTH | Mohakatino_River | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | MohakatinoRiver | 100.000000 |
| 10 | SOUTH | CapeEgmont | 0.58 | 0.78 | 0.84 | 1.10 | 0.0058 | 0.0078 | 0.0084 | 0.0110 | CapeEgmont | 100.000000 |
| 9 | NORTH | Waitara | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Waitara | 100.000000 |
| 8 | NORTH | UrenuiRiver_North | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | UrenuiRiverNorthCliffs | 100.000000 |
| 6 | NORTH | PariokariwaPoint | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | PariokariwaPointCliffs | 100.000000 |
| 5 | NORTH | Onaero | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | OnaeroCliff | 100.000000 |
| 4 | NORTH | Oakura_South | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | OakuraSouth | 100.000000 |
| 3 | NORTH | Oakura | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Oakura | 100.000000 |
| 2 | NORTH | NewPlymouth_Airport | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | NewPlymouthAirport | 100.000000 |
| 1 | NORTH | NewPlymouth | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | NewPlymouthCliffs | 100.000000 |
| 22 | SOUTH | WainuiBeach | 0.57 | 0.78 | 0.83 | 1.09 | 0.0057 | 0.0078 | 0.0083 | 0.0109 | WainuiBeach | 100.000000 |
| 23 | SOUTH | WaipipiBeach | 0.57 | 0.78 | 0.83 | 1.09 | 0.0057 | 0.0078 | 0.0083 | 0.0109 | WaipipiBeachCliffs | 100.000000 |
site = df.match.sample(1).iloc[0]
site
'UrenuiRiverNorthCliffs'
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
if site == "ManaBayCliffs":
print("Flipping")
for k, v in transect_metadata.items():
transect_metadata[k]["Azimuth"] = v["Azimuth"] + 180
results = predict(gdf, transect_metadata)
results
| TransectID | BaselineID | group | Year | ocean_point | linear_model_point | sqrt_model_point | BH_model_point | Sunamura_model_point | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 5 | 1 | 1 | 2100 | POINT (1720633.1100425983 5684136.280603779) | POINT (1720854.3339775638 5683672.05802439) | POINT (1720868.1061666948 5683643.158072117) | POINT (1720869.6189754275 5683639.983551366) | POINT (1720854.3339775638 5683672.05802439) |
| 1 | 6 | 1 | 1 | 2100 | POINT (1720642.0801050563 5684139.860357197) | POINT (1720862.6898028378 5683676.926710536) | POINT (1720873.4603526713 5683654.325483295) | POINT (1720877.9748007017 5683644.852237511) | POINT (1720862.6898028378 5683676.926710536) |
| 2 | 7 | 1 | 1 | 2100 | POINT (1720650.575075268 5684144.8879364645) | POINT (1720871.3119035014 5683681.687515657) | POINT (1720881.2773631075 5683660.775712584) | POINT (1720886.5969013653 5683649.613042634) | POINT (1720871.3119035014 5683681.687515657) |
| 3 | 8 | 1 | 1 | 2100 | POINT (1720659.1581115604 5684150.087639844) | POINT (1720879.2862932794 5683688.164420146) | POINT (1720888.6494842013 5683668.51643476) | POINT (1720894.571291143 5683656.089947121) | POINT (1720879.2862932794 5683688.164420146) |
| 4 | 9 | 1 | 1 | 2100 | POINT (1720667.0028565435 5684155.18472246) | POINT (1720886.9824249977 5683693.57335705) | POINT (1720897.615215984 5683671.261206882) | POINT (1720902.2674228614 5683661.498884027) | POINT (1720886.9824249977 5683693.57335705) |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 314 | 414 | 25 | 31 | 2100 | POINT (1726658.0317596905 5689964.909653514) | POINT (1727097.8887274514 5689705.400302266) | POINT (1727116.2523344012 5689694.5660355035) | POINT (1727128.490082377 5689687.345938357) | POINT (1727097.8887274514 5689705.400302266) |
| 315 | 415 | 25 | 31 | 2100 | POINT (1726662.3798315164 5689973.425984194) | POINT (1727101.3898836626 5689714.932497662) | POINT (1727116.3902041595 5689706.100159711) | POINT (1727132.0069823388 5689696.904845344) | POINT (1727101.3898836626 5689714.932497662) |
| 316 | 416 | 25 | 31 | 2100 | POINT (1726666.1429697142 5689982.1815043455) | POINT (1727103.4452041262 5689725.309927042) | POINT (1727114.616595727 5689718.747845362) | POINT (1727134.0811645247 5689707.314346753) | POINT (1727103.4452041262 5689725.309927042) |
| 317 | 417 | 25 | 31 | 2100 | POINT (1726669.315968925 5689991.279733421) | POINT (1727107.1794726513 5689734.6948276665) | POINT (1727119.1615009448 5689727.67344412) | POINT (1727137.834261176 5689716.731339142) | POINT (1727107.1794726513 5689734.6948276665) |
| 318 | 418 | 25 | 31 | 2100 | POINT (1726672.9711513557 5690000.091062417) | POINT (1727113.035012376 5689742.835461913) | POINT (1727129.9653122795 5689732.938229137) | POINT (1727143.70859541 5689724.904084851) | POINT (1727113.035012376 5689742.835461913) |
319 rows × 9 columns
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
poly = prediction_results_to_polygon(results)
poly.explore(tiles="Esri.WorldImagery")
site = df.match.sample(1).iloc[0]
site = "CapeEgmont"
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
m = gpd.read_file(f"Shapefiles/{site}_TransectLines.shp").explore("TransectID", tiles="Esri.WorldImagery", max_zoom=22)
results = predict(gdf, transect_metadata)
results["geometry"] = results["linear_model_point"]
pd.concat([gdf, results])[["Year", "geometry"]].explore("Year", m=m)
gpd.GeoSeries(results.ocean_point, crs=2193).explore(m=m, color="blue")
#gpd.GeoSeries(poly, crs=2193).explore(color="pink", m=m)
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
polygon = prediction_results_to_polygon(results)
polygon.explore(m=m)
m
CapeEgmont
run_all_parallel()
0%| | 0/24 [00:00<?, ?it/s]
Flipping
def read_file(f):
df = gpd.read_file(f)
df["filename"] = f
return df
samples = pd.concat(read_file(f) for f in glob("Projected_Shoreline_Polygons/*.shp"))
samples.explore("filename", tiles="Esri.WorldImagery", style_kwds={"fill": False})